module gw_drag

  !--------------------------------------------------------------------------
  ! CAM and WACCM gravity wave parameterizations were merged by Sean Patrick
  ! Santos in Summer 2013, and at the same time, gw_drag was split into
  ! various modules. This is the CAM interface and driver module. The below
  ! notes are for the old CAM and WACCM versions of gw_drag.
  !--------------------------------------------------------------------------
  ! This file came from wa17 and was modified by Fabrizio: 07-02-2004
  ! Standard gw_drag with modification (6) of latitude profile of gw spectrum
  !--------------------------------------------------------------------------
  ! Purpose:
  !
  ! Module to compute the forcing due to parameterized gravity waves. Both an
  ! orographic and an internal source spectrum are considered.
  !
  ! Author: Byron Boville
  !--------------------------------------------------------------------------

  use shr_kind_mod  , only: r8 => shr_kind_r8, cl => shr_kind_cl
  use shr_log_mod   , only: errMsg => shr_log_errMsg
  use shr_assert_mod, only: shr_assert
  use ppgrid        , only: pcols, pver, begchunk, endchunk
  use constituents  , only: pcnst
  use physics_types , only: physics_state, physics_ptend, physics_ptend_init
  use spmd_utils    , only: masterproc
  use cam_history   , only: outfld
  use cam_logfile   , only: iulog
  use cam_abortutils, only: endrun
  use ref_pres      , only: do_molec_diff, nbot_molec, press_lim_idx
  use physconst     , only: cpair

  ! These are the actual switches for different gravity wave sources.
  use phys_control  , only: use_gw_oro, use_gw_front, use_gw_front_igw, &
                            use_gw_convect_dp, use_gw_convect_sh,       &
                            use_simple_phys
  use gw_common     , only: GWBand
  use gw_convect    , only: BeresSourceDesc
  use gw_front      , only: CMSourceDesc

  implicit none

  private

  save

  public gw_drag_readnl           ! Read namelist
  public gw_init                  ! Initialization
  public gw_tend                  ! interface to actual parameterization

  real(r8), parameter :: unset_r8 = huge(1.0_r8)

  ! A mid-scale "band" with only stationary waves (l = 0).
  type(GWBand) band_oro
  ! Medium scale waves.
  type(GWBand) band_mid
  ! Long scale waves for IGWs.
  type(GWBand) band_long

  ! Top level for gravity waves.
  integer, parameter :: ktop = 1
  ! Bottom level for frontal waves.
  integer kbot_front

  ! Factor for SH orographic waves.
  real(r8) :: gw_oro_south_fac = 1.0_r8

  ! Frontogenesis function critical threshold.
  real(r8) :: frontgfc = unset_r8

  ! Tendency efficiencies.

  ! Ridge scheme.
  logical  :: use_gw_rdg_beta      = .false.
  integer  :: n_rdg_beta           = 1
  real(r8) :: effgw_rdg_beta       = unset_r8
  real(r8) :: effgw_rdg_beta_max   = unset_r8
  real(r8) :: rdg_beta_cd_llb      = unset_r8  ! Low-level obstacle drag coefficient Ridge scheme.
  logical  :: trpd_leewv_rdg_beta  = .false.

  logical  :: use_gw_rdg_gamma     = .false.
  integer  :: n_rdg_gamma          = -1
  real(r8) :: effgw_rdg_gamma      = unset_r8
  real(r8) :: effgw_rdg_gamma_max  = unset_r8
  real(r8) :: rdg_gamma_cd_llb     = unset_r8
  logical  :: trpd_leewv_rdg_gamma = .false.
  character(cl) :: bnd_rdggm       = 'bnd_rdggm' ! full pathname for meso-Gamma ridge dataset

  ! Orography.
  real(r8) :: effgw_oro = unset_r8
  ! C&M scheme.
  real(r8) :: effgw_cm = unset_r8
  ! C&M scheme (inertial waves).
  real(r8) :: effgw_cm_igw = unset_r8
  ! Beres (deep convection).
  real(r8) :: effgw_beres_dp = unset_r8
  ! Beres (shallow convection).
  real(r8) :: effgw_beres_sh = unset_r8

  ! Horzontal wavelengths [m].
  real(r8), parameter :: wavelength_mid  = 1.e5_r8
  real(r8), parameter :: wavelength_long = 3.e5_r8

  ! Background stress source strengths.
  real(r8) :: taubgnd     = unset_r8
  real(r8) :: taubgnd_igw = unset_r8

  ! Whether or not to use a polar taper for frontally generated waves.
  logical :: gw_polar_taper = .false.

  ! Whether or not to enforce an upper boundary condition of tau = 0.
  ! (Like many variables, this is only here to hold the value between
  ! the readnl phase and the init phase of the CAM physics; only gw_common
  ! should actually use it.)
  logical :: tau_0_ubc = .false.

  ! Whether or not to limit tau *before* applying any efficiency factors.
  logical :: gw_limit_tau_without_eff = .false.

  ! Whether or not to apply tendency max
  logical :: gw_apply_tndmax = .true.

  ! Files to read Beres source spectra from.
  character(256) :: gw_drag_file = ""
  character(256) :: gw_drag_file_sh = ""

  ! Beres settings and table.
  type(BeresSourceDesc) beres_dp_desc
  type(BeresSourceDesc) beres_sh_desc

  ! Width of gaussian used to create frontogenesis tau profile [m/s].
  real(r8), parameter :: front_gaussian_width = 30.0_r8

  ! Frontogenesis wave settings.
  type(CMSourceDesc) cm_desc
  type(CMSourceDesc) cm_igw_desc

  ! Indices into pbuf
  integer :: kvt_idx      = -1
  integer :: ttend_dp_idx = -1
  integer :: ttend_sh_idx = -1
  integer :: frontgf_idx  = -1
  integer :: frontga_idx  = -1
  integer :: sgh_idx      = -1

  ! anisotropic ridge fields
  integer, parameter :: prdg = 16

  real(r8), allocatable, dimension(:,:), target :: rdg_gbxar

  ! Meso Beta
  real(r8), allocatable, dimension(:,:,:), target :: &
    rdg_hwdth,  &
    rdg_clngt,  &
    rdg_mxdis,  &
    rdg_anixy,  &
    rdg_angll

  ! Meso Gamma
  real(r8), allocatable, dimension(:,:,:), target :: &
    rdg_hwdthg, &
    rdg_clngtg, &
    rdg_mxdisg, &
    rdg_anixyg, &
    rdg_angllg

  ! Water constituent indices for budget
  integer :: ixcldliq = -1
  integer :: ixcldice = -1

  ! Prefixes for history field names
  character(1), parameter :: cm_pf = " "
  character(1), parameter :: cm_igw_pf = "I"
  character(1), parameter :: beres_dp_pf = "B"
  character(1), parameter :: beres_sh_pf = "S"

  ! namelist
  logical  :: history_amwg                    ! Output the variables used by the AMWG diag package
  logical  :: gw_lndscl_sgh         = .true.  ! Scale SGH by land fraction
  real(r8) :: gw_prndl              = 0.25_r8
  real(r8) :: gw_qbo_hdepth_scaling = 1.0_r8  ! Heating depth scaling factor

  logical :: gw_top_taper = .false.
  real(r8), pointer :: vramp(:) => null()

contains

  subroutine gw_drag_readnl(nlfile)

    use namelist_utils, only: find_group_name
    use units,          only: getunit, freeunit
    use spmd_utils,     only: mpicom, mstrid=>masterprocid, mpi_real8, &
                              mpi_character, mpi_logical, mpi_integer
    use gw_rdg,         only: gw_rdg_readnl

    character(*), intent(in) :: nlfile

    integer unitn, ierr
    character(*), parameter :: sub = 'gw_drag_readnl'

    ! Maximum wave number and width of spectrum bins.
    integer  :: pgwv       = -1
    real(r8) :: gw_dc      = unset_r8
    integer  :: pgwv_long  = -1
    real(r8) :: gw_dc_long = unset_r8

    ! fcrit2 for the mid-scale waves has been made a namelist variable to
    ! facilitate backwards compatibility with the CAM3 version of this
    ! parameterization.  In CAM3, fcrit2=0.5.
    real(r8) :: fcrit2 = unset_r8   ! critical froude number squared

    namelist /gw_drag_nl/ pgwv, gw_dc, pgwv_long, gw_dc_long, tau_0_ubc, &
      effgw_beres_dp, effgw_beres_sh, effgw_cm, effgw_cm_igw, effgw_oro, &
      fcrit2, frontgfc, gw_drag_file, gw_drag_file_sh, taubgnd, &
      taubgnd_igw, gw_polar_taper, &
      use_gw_rdg_beta, n_rdg_beta, effgw_rdg_beta, effgw_rdg_beta_max, &
      rdg_beta_cd_llb, trpd_leewv_rdg_beta, &
      use_gw_rdg_gamma, n_rdg_gamma, effgw_rdg_gamma, effgw_rdg_gamma_max, &
      rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, bnd_rdggm, &
      gw_oro_south_fac, gw_limit_tau_without_eff, &
      gw_lndscl_sgh, gw_prndl, gw_apply_tndmax, gw_qbo_hdepth_scaling, &
      gw_top_taper

    if (use_simple_phys) return

    if (masterproc) then
      unitn = getunit()
      open(unitn, file=trim(nlfile), status='old')
      call find_group_name(unitn, 'gw_drag_nl', status=ierr)
      if (ierr == 0) then
        read(unitn, gw_drag_nl, iostat=ierr)
        if (ierr /= 0) then
          call endrun(sub // ':: ERROR reading namelist')
        end if
      end if
      close(unitn)
      call freeunit(unitn)
    end if

    call mpi_bcast(pgwv, 1, mpi_integer, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: pgwv")
    call mpi_bcast(gw_dc, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_dc")
    call mpi_bcast(pgwv_long, 1, mpi_integer, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: pgwv_long")
    call mpi_bcast(gw_dc_long, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_dc_long")
    call mpi_bcast(tau_0_ubc, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: tau_0_ubc")
    call mpi_bcast(effgw_beres_dp, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_beres_dp")
    call mpi_bcast(effgw_beres_sh, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_beres_sh")
    call mpi_bcast(effgw_cm, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_cm")
    call mpi_bcast(effgw_cm_igw, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_cm_igw")
    call mpi_bcast(effgw_oro, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_oro")

    call mpi_bcast(use_gw_rdg_beta, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_gw_rdg_beta")
    call mpi_bcast(n_rdg_beta, 1, mpi_integer, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: n_rdg_beta")
    call mpi_bcast(effgw_rdg_beta, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_beta")
    call mpi_bcast(effgw_rdg_beta_max, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_beta_max")
    call mpi_bcast(rdg_beta_cd_llb, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rdg_beta_cd_llb")
    call mpi_bcast(trpd_leewv_rdg_beta, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: trpd_leewv_rdg_beta")

    call mpi_bcast(use_gw_rdg_gamma, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_gw_rdg_gamma")
    call mpi_bcast(n_rdg_gamma, 1, mpi_integer, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: n_rdg_gamma")
    call mpi_bcast(effgw_rdg_gamma, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_gamma")
    call mpi_bcast(effgw_rdg_gamma_max, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_gamma_max")
    call mpi_bcast(rdg_gamma_cd_llb, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rdg_gamma_cd_llb")
    call mpi_bcast(trpd_leewv_rdg_gamma, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: trpd_leewv_rdg_gamma")
    call mpi_bcast(bnd_rdggm, len(bnd_rdggm), mpi_character, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: bnd_rdggm")

    call mpi_bcast(gw_oro_south_fac, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_oro_south_fac")
    call mpi_bcast(fcrit2, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: fcrit2")
    call mpi_bcast(frontgfc, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frontgfc")
    call mpi_bcast(taubgnd, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: taubgnd")
    call mpi_bcast(taubgnd_igw, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: taubgnd_igw")

    call mpi_bcast(gw_polar_taper, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_polar_taper")
    call mpi_bcast(gw_limit_tau_without_eff, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_limit_tau_without_eff")
    call mpi_bcast(gw_apply_tndmax, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_apply_tndmax")

    call mpi_bcast(gw_top_taper, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_top_taper")

    call mpi_bcast(gw_lndscl_sgh, 1, mpi_logical, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_lndscl_sgh")
    call mpi_bcast(gw_prndl, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_prndl")
    call mpi_bcast(gw_qbo_hdepth_scaling, 1, mpi_real8, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_qbo_hdepth_scaling")

    call mpi_bcast(gw_drag_file, len(gw_drag_file), mpi_character, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_drag_file")
    call mpi_bcast(gw_drag_file_sh, len(gw_drag_file_sh), mpi_character, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_drag_file_sh")


    ! Check if fcrit2 was set.
    call shr_assert(fcrit2 /= unset_r8, 'gw_drag_readnl: fcrit2 must be set via the namelist.' // errMsg(__FILE__, __LINE__))

    ! Check if pgwv was set.
    call shr_assert(pgwv >= 0, 'gw_drag_readnl: pgwv must be set via the namelist and non-negative.' // errMsg(__FILE__, __LINE__))

    ! Check if gw_dc was set.
    call shr_assert(gw_dc /= unset_r8, 'gw_drag_readnl: gw_dc must be set via the namelist.' // errMsg(__FILE__, __LINE__))

    band_oro = GWBand(0, gw_dc, fcrit2, wavelength_mid)
    band_mid = GWBand(pgwv, gw_dc, 1.0_r8, wavelength_mid)
    band_long = GWBand(pgwv_long, gw_dc_long, 1.0_r8, wavelength_long)

    if (use_gw_rdg_gamma .or. use_gw_rdg_beta) call gw_rdg_readnl(nlfile)

  end subroutine gw_drag_readnl

  subroutine gw_init()

    !-----------------------------------------------------------------------
    ! Time independent initialization for multiple gravity wave
    ! parameterization.
    !-----------------------------------------------------------------------

    use cam_history,      only: addfld, add_default, horiz_only
    use cam_history,      only: register_vector_field
    use interpolate_data, only: lininterp
    use phys_control,     only: phys_getopts
    use physics_buffer,   only: pbuf_get_index
    use constituents,     only: cnst_get_ind

    use cam_initfiles,    only: topo_file_get_id

    ! temporary for restart with ridge scheme
    use cam_initfiles,    only: bnd_topo

    use cam_pio_utils,    only: cam_pio_openfile
    use cam_grid_support, only: cam_grid_check, cam_grid_id
    use cam_grid_support, only: cam_grid_get_dim_names
    use pio,              only: file_desc_t, pio_nowrite, pio_closefile
    use ncdio_atm,        only: infld
    use ioFileMod,        only: getfil

    use ref_pres,         only: pref_edge, pref_mid
    use physconst,        only: gravit, rair, rearth

    use gw_common,        only: gw_common_init
    use gw_front,         only: gaussian_cm_desc

    integer i, l, k
    character(1) cn

    ! Index for levels at specific pressures.
    integer kfront

    ! output tendencies and state variables for CAM4 temperature,
    ! water vapor, cloud ice and cloud liquid budgets.
    logical history_budget
    ! output history file number for budget fields
    integer history_budget_histfile_num
    ! output variables of interest in WACCM runs
    logical history_waccm

    ! Interpolated Newtonian cooling coefficients.
    real(r8) alpha(pver+1)

    ! Levels of pre-calculated Newtonian cooling (1/day).
    ! The following profile is digitized from:
    ! Wehrbein and Leovy (JAS, 39, 1532-1544, 1982) figure 5
    integer, parameter :: nalph = 71
    real(r8) :: alpha0(nalph) = [ &
      0.1_r8,         0.1_r8,         0.1_r8,         0.1_r8,         &
      0.1_r8,         0.1_r8,         0.1_r8,         0.1_r8,         &
      0.1_r8,         0.1_r8,         0.10133333_r8,  0.104_r8,       &
      0.108_r8,       0.112_r8,       0.116_r8,       0.12066667_r8,  &
      0.126_r8,       0.132_r8,       0.138_r8,       0.144_r8,       &
      0.15133333_r8,  0.16_r8,        0.17_r8,        0.18_r8,        &
      0.19_r8,        0.19933333_r8,  0.208_r8,       0.216_r8,       &
      0.224_r8,       0.232_r8,       0.23466667_r8,  0.232_r8,       &
      0.224_r8,       0.216_r8,       0.208_r8,       0.20133333_r8,  &
      0.196_r8,       0.192_r8,       0.188_r8,       0.184_r8,       &
      0.18266667_r8,  0.184_r8,       0.188_r8,       0.192_r8,       &
      0.196_r8,       0.19333333_r8,  0.184_r8,       0.168_r8,       &
      0.152_r8,       0.136_r8,       0.12133333_r8,  0.108_r8,       &
      0.096_r8,       0.084_r8,       0.072_r8,       0.061_r8,       &
      0.051_r8,       0.042_r8,       0.033_r8,       0.024_r8,       &
      0.017666667_r8, 0.014_r8,       0.013_r8,       0.012_r8,       &
      0.011_r8,       0.010333333_r8, 0.01_r8,        0.01_r8,        &
      0.01_r8,        0.01_r8,        0.01_r8                         &
    ]

    ! Pressure levels that were used to calculate alpha0 (hPa).
    real(r8) :: palph(nalph) = [ &
      2.06115E-06_r8, 2.74280E-06_r8, 3.64988E-06_r8, 4.85694E-06_r8, &
      6.46319E-06_r8, 8.60065E-06_r8, 1.14450E-05_r8, 1.52300E-05_r8, &
      2.02667E-05_r8, 2.69692E-05_r8, 3.58882E-05_r8, 4.77568E-05_r8, &
      6.35507E-05_r8, 8.45676E-05_r8, 0.000112535_r8, 0.000149752_r8, &
      0.000199277_r8, 0.000265180_r8, 0.000352878_r8, 0.000469579_r8, &
      0.000624875_r8, 0.000831529_r8, 0.00110653_r8,  0.00147247_r8,  &
      0.00195943_r8,  0.00260744_r8,  0.00346975_r8,  0.00461724_r8,  &
      0.00614421_r8,  0.00817618_r8,  0.0108801_r8,   0.0144783_r8,   &
      0.0192665_r8,   0.0256382_r8,   0.0341170_r8,   0.0453999_r8,   &
      0.0604142_r8,   0.0803939_r8,   0.106981_r8,    0.142361_r8,    &
      0.189442_r8,    0.252093_r8,    0.335463_r8,    0.446404_r8,    &
      0.594036_r8,    0.790490_r8,    1.05192_r8,     1.39980_r8,     &
      1.86273_r8,     2.47875_r8,     3.29851_r8,     4.38936_r8,     &
      5.84098_r8,     7.77266_r8,     10.3432_r8,     13.7638_r8,     &
      18.3156_r8,     24.3728_r8,     32.4332_r8,     43.1593_r8,     &
      57.4326_r8,     76.4263_r8,     101.701_r8,     135.335_r8,     &
      180.092_r8,     239.651_r8,     318.907_r8,     424.373_r8,     &
      564.718_r8,     751.477_r8,     1000.0_r8                        &
    ]

    ! Read data from file
    type(file_desc_t), pointer :: fh_topo
    type(file_desc_t)  fh_rdggm
    integer            grid_id
    character(8)       dim1name, dim2name
    logical            found
    character(256)     bnd_rdggm_loc   ! File path of topo file on local disk

    ! Allow reporting of error messages.
    character(128)     errstring
    character(*), parameter :: sub = 'gw_init'

    ! temporary workaround for restart w/ ridge scheme
    character(256)     bnd_topo_loc   ! File path of topo file on local disk

    integer botndx,topndx

    if (do_molec_diff) then
      kvt_idx = pbuf_get_index('kvt')
    end if

    if (masterproc) then
      write(iulog, *) ' '
      write(iulog, *) 'GW_DRAG: band_mid%ngwv = ', band_mid%ngwv
      do l = -band_mid%ngwv, band_mid%ngwv
          write(iulog,'(A,I0,A,F7.2)') ' GW_DRAG: band_mid%cref(', l, ') = ', band_mid%cref(l)
      end do
      write(iulog, *) 'GW_DRAG: band_mid%kwv = ', band_mid%kwv
      write(iulog, *) 'GW_DRAG: band_mid%fcrit2 = ', band_mid%fcrit2
      write(iulog, *) ' '
      write(iulog, *) 'GW_DRAG: band_long%ngwv = ', band_long%ngwv
      do l = -band_long%ngwv, band_long%ngwv
          write(iulog,'(A,I2,A,F7.2)') ' GW_DRAG: band_long%cref(', l, ') = ', band_long%cref(l)
      end do
      write(iulog, *) 'GW_DRAG: band_long%kwv = ', band_long%kwv
      write(iulog, *) 'GW_DRAG: band_long%fcrit2 = ', band_long%fcrit2
      write(iulog, *) ' '
    end if

    ! pre-calculated newtonian damping:
    !     * convert to 1/s
    !     * ensure it is not smaller than 1e-6
    !     * convert palph from hpa to pa

    do k = 1, nalph
      alpha0(k) = alpha0(k) / 86400.0_r8
      alpha0(k) = max(alpha0(k), 1.0e-6_r8)
      palph (k) = palph(k) * 1.0e2_r8
    end do

    ! Interpolate to current vertical grid and obtain alpha
    call lininterp(alpha0  ,palph, nalph , alpha  , pref_edge , pver+1)
    if (masterproc) then
      write(iulog, *) 'gw_init: newtonian damping (1/day):'
      write(iulog, '(A4, A12, A10)') ' k  ','  pref_edge      ', '  alpha   '
      do k = 1, pver + 1
        write(iulog, '(I4, 1e12.5, 1f10.2)') k, pref_edge(k), alpha(k) * 86400
      end do
    end if

    if (masterproc) then
      write(iulog, *) 'KTOP        =', ktop
    end if

    ! Used to decide whether temperature tendencies should be output.
    call phys_getopts(                                                &
        history_budget_out              =history_budget             , &
        history_budget_histfile_num_out =history_budget_histfile_num, &
        history_waccm_out               =history_waccm              , &
        history_amwg_out                =history_amwg               )

    ! Initialize subordinate modules.
    call gw_common_init(pver, tau_0_ubc, ktop, gravit, rair, alpha, gw_prndl, &
                        gw_qbo_hdepth_scaling, errstring)
    call shr_assert(trim(errstring) == '', 'gw_common_init: ' // errstring // &
                    errMsg(__FILE__, __LINE__))

    if (use_gw_oro) then
      if (effgw_oro == unset_r8) then
        call endrun('gw_drag_init: Orographic gravity waves enabled, but effgw_oro was not set.')
      end if
    end if

    if (use_gw_oro .or. use_gw_rdg_beta .or. use_gw_rdg_gamma) then
      sgh_idx = pbuf_get_index('SGH')
      ! Declare history variables for orographic term
      call addfld ('TAUAORO'   , ['ilev'  ], 'I', 'N/m2', 'Total stress from original OGW scheme')
      call addfld ('TTGWORO'   , ['lev'   ], 'A', 'K/s' , 'T tendency - orographic gravity wave drag')
      call addfld ('TTGWSDFORO', ['lev'   ], 'A', 'K/s' , 'T tendency - orographic gravity wave, diffusion.')
      call addfld ('TTGWSKEORO', ['lev'   ], 'A', 'K/s' , 'T tendency - orographic gravity wave, breaking KE.')
      call addfld ('UTGWORO'   , ['lev'   ], 'A', 'm/s2', 'U tendency - orographic gravity wave drag')
      call addfld ('VTGWORO'   , ['lev'   ], 'A', 'm/s2', 'V tendency - orographic gravity wave drag')
      call register_vector_field('UTGWORO', 'VTGWORO')
      call addfld ('TAUGWX'    , horiz_only, 'A', 'N/m2', 'Zonal gravity wave surface stress')
      call addfld ('TAUGWY'    , horiz_only, 'A', 'N/m2', 'Meridional gravity wave surface stress')
      call register_vector_field('TAUGWX', 'TAUGWY')

      if (history_amwg) then
        call add_default('TAUGWX  ', 1, ' ')
        call add_default('TAUGWY  ', 1, ' ')
      end if

      if (history_waccm) then
        call add_default('UTGWORO ', 1, ' ')
        call add_default('VTGWORO ', 1, ' ')
        call add_default('TAUGWX  ', 1, ' ')
        call add_default('TAUGWY  ', 1, ' ')
      end if
    end if

    if (use_gw_rdg_beta .or. use_gw_rdg_gamma) then
      grid_id = cam_grid_id('physgrid')
      if (.not. cam_grid_check(grid_id)) then
        call endrun(sub // ': ERROR: no "physgrid" grid')
      end if
      call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
    end if

    if (use_gw_rdg_beta) then
      if (effgw_rdg_beta == unset_r8) then
        call endrun(sub // ': ERROR: Anisotropic OGW enabled, but effgw_rdg_beta was not set.')
      end if

      fh_topo => topo_file_get_id()
      bnd_topo_loc = ' '
      if (.not. associated(fh_topo)) then
        ! Try to open topo file here.  This workaround will not be needed
        ! once the refactored initialization sequence is on trunk.
        allocate(fh_topo)
        ! Error exit is from getfil if file not found.
        call getfil(bnd_topo, bnd_topo_loc)
        call cam_pio_openfile(fh_topo, bnd_topo_loc, PIO_NOWRITE)
      end if

      ! Get beta ridge data
      allocate( &
        rdg_gbxar(pcols,     begchunk:endchunk), &
        rdg_hwdth(pcols,prdg,begchunk:endchunk), &
        rdg_clngt(pcols,prdg,begchunk:endchunk), &
        rdg_mxdis(pcols,prdg,begchunk:endchunk), &
        rdg_anixy(pcols,prdg,begchunk:endchunk), &
        rdg_angll(pcols,prdg,begchunk:endchunk))

      call infld('GBXAR', fh_topo, dim1name, dim2name, 1, pcols, &
                 begchunk, endchunk, rdg_gbxar, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: GBXAR not found on topo file')
      rdg_gbxar = rdg_gbxar * (rearth/1000.0_r8)*(rearth/1000.0_r8) ! transform to km^2

      call infld('HWDTH', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
                 1, prdg, begchunk, endchunk, rdg_hwdth, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: HWDTH not found on topo file')

      call infld('CLNGT', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
                 1, prdg, begchunk, endchunk, rdg_clngt, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: CLNGT not found on topo file')

      call infld('MXDIS', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
                 1, prdg, begchunk, endchunk, rdg_mxdis, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: MXDIS not found on topo file')

      call infld('ANIXY', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
                 1, prdg, begchunk, endchunk, rdg_anixy, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: ANIXY not found on topo file')

      call infld('ANGLL', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
                 1, prdg, begchunk, endchunk, rdg_angll, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: ANGLL not found on topo file')

      ! close topo file only if it was opened here
      if (len_trim(bnd_topo_loc) > 0) then
        call pio_closefile(fh_topo)
      end if

      call addfld('UEGW'        , ['lev' ], 'A', '1/s' , 'Zonal wind profile-entry to GW')
      call addfld('VEGW'        , ['lev' ], 'A', '1/s' , 'Merdional wind profile-entry to GW')
      call register_vector_field('UEGW','VEGW')
      call addfld('TEGW'        , ['lev' ], 'A', 'K'   , 'Temperature profile-entry to GW')
      call addfld('TAU1RDGBETAM', ['ilev'], 'I', 'N/m2', 'Ridge based momentum flux profile')
      call addfld('UBM1BETA'    , ['lev' ], 'A', '1/s' , 'On-ridge wind profile')
      call addfld('UBT1RDGBETA' , ['lev' ], 'I', 'm/s' , 'On-ridge wind tendency from ridge 1')

      do i = 1, 6
        write(cn, '(I1)') i
        call addfld('TAU' // cn // 'RDGBETAY', ['ilev'], 'I', 'N/m2', 'Ridge based momentum flux profile')
        call addfld('TAU' // cn // 'RDGBETAX', ['ilev'], 'I', 'N/m2', 'Ridge based momentum flux profile')
        call register_vector_field('TAU' // cn // 'RDGBETAX', 'TAU' // cn // 'RDGBETAY')
        call addfld('UT'  // cn // 'RDGBETA' , ['lev' ], 'I', 'm/s' , 'U wind tendency from ridge ' // cn)
        call addfld('VT'  // cn // 'RDGBETA' , ['lev' ], 'I', 'm/s' , 'V wind tendency from ridge ' // cn)
        call register_vector_field('UT' // cn // 'RDGBETA', 'VT' // cn // 'RDGBETA')
      end do

      call addfld('TAUARDGBETAY', ['ilev'], 'I', 'N/m2', 'Ridge based momentum flux profile')
      call addfld('TAUARDGBETAX', ['ilev'], 'I', 'N/m2', 'Ridge based momentum flux profile')
      call register_vector_field('TAUARDGBETAX','TAUARDGBETAY')

      if (history_waccm) then
        call add_default('TAUARDGBETAX', 1, ' ')
        call add_default('TAUARDGBETAY', 1, ' ')
      end if
    end if

    if (use_gw_rdg_gamma) then
      if (effgw_rdg_gamma == unset_r8) then
        call endrun(sub // ': ERROR: Anisotropic OGW enabled, but effgw_rdg_gamma was not set.')
      end if

      call getfil(bnd_rdggm, bnd_rdggm_loc, iflag=1, lexist=found)
      if (found) then
        call cam_pio_openfile(fh_rdggm, bnd_rdggm_loc, PIO_NOWRITE)
      else
        call endrun(sub // ': ERROR: file for gamma ridges not found: bnd_rdggm=' // trim(bnd_rdggm))
      end if

      if (.not. allocated(rdg_gbxar)) then
        allocate(rdg_gbxar(pcols,begchunk:endchunk))
        call infld('GBXAR', fh_rdggm, dim1name, dim2name, 1, pcols, &
                   begchunk, endchunk, rdg_gbxar, found, gridname='physgrid')
        if (.not. found) call endrun(sub // ': ERROR: GBXAR not found on bnd_rdggm')
        rdg_gbxar = rdg_gbxar * (rearth / 1000.0_r8)**2 ! transform to km^2
      end if

      ! Get meso-gamma ridge data
      allocate( &
        rdg_hwdthg(pcols,prdg,begchunk:endchunk), &
        rdg_clngtg(pcols,prdg,begchunk:endchunk), &
        rdg_mxdisg(pcols,prdg,begchunk:endchunk), &
        rdg_anixyg(pcols,prdg,begchunk:endchunk), &
        rdg_angllg(pcols,prdg,begchunk:endchunk))

      call infld('HWDTH', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
                  1, prdg, begchunk, endchunk, rdg_hwdthg, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: HWDTH not found on bnd_rdggm')

      call infld('CLNGT', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
                  1, prdg, begchunk, endchunk, rdg_clngtg, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: CLNGT not found on bnd_rdggm')

      call infld('MXDIS', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
                  1, prdg, begchunk, endchunk, rdg_mxdisg, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: MXDIS not found on bnd_rdggm')

      call infld('ANIXY', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
                  1, prdg, begchunk, endchunk, rdg_anixyg, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: ANIXY not found on bnd_rdggm')

      call infld('ANGLL', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
                  1, prdg, begchunk, endchunk, rdg_angllg, found, gridname='physgrid')
      if (.not. found) call endrun(sub//': ERROR: ANGLL not found on bnd_rdggm')

      call pio_closefile(fh_rdggm)

      call addfld('TAU1RDGGAMMAM', ['ilev'], 'I', 'N/m2', 'Ridge based momentum flux profile')
      call addfld('UBM1GAMMA'    , ['lev' ], 'A', '1/s' , 'On-ridge wind profile')
      call addfld('UBT1RDGGAMMA' , ['lev' ], 'I', 'm/s' , 'On-ridge wind tendency from ridge 1')

      do i = 1, 6
        write(cn, '(i1)') i
        call addfld('TAU' // cn // 'RDGGAMMAY', ['ilev'], 'I', 'N/m2', 'Ridge based momentum flux profile')
        call addfld('TAU' // cn // 'RDGGAMMAX', ['ilev'], 'I', 'N/m2', 'Ridge based momentum flux profile')
        call addfld('UT'  // cn // 'RDGGAMMA' , ['lev' ], 'I', 'm/s' , 'U wind tendency from ridge ' // cn)
        call addfld('VT'  // cn // 'RDGGAMMA' , ['lev' ], 'I', 'm/s' , 'V wind tendency from ridge ' // cn)
        call register_vector_field('UT' // cn // 'RDGGAMMA', 'VT' // cn // 'RDGGAMMA')
      end do

      call addfld('TAUARDGGAMMAY', [ 'ilev' ], 'I', 'N/m2', 'Ridge based momentum flux profile')
      call addfld('TAUARDGGAMMAX', [ 'ilev' ], 'I', 'N/m2', 'Ridge based momentum flux profile')
      call register_vector_field('TAUARDGGAMMAX','TAUARDGGAMMAY')
      call addfld('TAURDGGMX', horiz_only, 'A', 'N/m2', 'Zonal gravity wave surface stress')
      call addfld('TAURDGGMY', horiz_only, 'A', 'N/m2', 'Meridional gravity wave surface stress')
      call register_vector_field('TAURDGGMX','TAURDGGMY')
      call addfld('UTRDGGM', ['lev'], 'I', 'm/s', 'U wind tendency from ridge 6')
      call addfld('VTRDGGM', ['lev'], 'I', 'm/s', 'V wind tendency from ridge 6')
      call register_vector_field('UTRDGGM','VTRDGGM')
    end if

    if (use_gw_front .or. use_gw_front_igw) then
      frontgf_idx = pbuf_get_index('FRONTGF')
      frontga_idx = pbuf_get_index('FRONTGA')
      call shr_assert(unset_r8 /= frontgfc, &
        'gw_drag_init: Frontogenesis enabled, but frontgfc was not set!' // errMsg(__FILE__, __LINE__))
      do k = 0, pver
        ! Check frontogenesis at 600 hPa.
        if (pref_edge(k+1) < 60000.0_r8) kfront = k + 1
      end do

      ! Source waves from 500 hPa.
      kbot_front = maxloc(pref_edge, 1, (pref_edge < 50000.0_r8)) - 1

      if (masterproc) then
        write(iulog, *) 'KFRONT      =', kfront
        write(iulog, *) 'KBOT_FRONT  =', kbot_front
        write(iulog, *) ' '
      end if

      call addfld('FRONTGF' , ['lev'], 'A', 'K^2/M^2/S', 'Frontogenesis function at gws src level')
      call addfld('FRONTGFA', ['lev'], 'A', 'K^2/M^2/S', 'Frontogenesis function at gws src level')

      if (history_waccm) then
        call add_default('FRONTGF' , 1, ' ')
        call add_default('FRONTGFA', 1, ' ')
      end if
    end if

    if (use_gw_front) then
      call shr_assert(all(unset_r8 /= [effgw_cm, taubgnd]), &
        'gw_drag_init: Frontogenesis mid-scale waves enabled, but not all required namelist variables were set!' // &
        errMsg(__FILE__, __LINE__))
      if (masterproc) then
        write(iulog, *) 'gw_init: gw spectrum taubgnd, effgw_cm = ', taubgnd, effgw_cm
        write(iulog, *) ' '
      end if
      cm_desc = gaussian_cm_desc(band_mid, kbot_front, kfront, frontgfc, taubgnd, front_gaussian_width)
      ! Output for gravity waves from frontogenesis.
      call gw_spec_addflds(prefix=cm_pf, scheme='C&M', band=band_mid, history_defaults=history_waccm)
    end if

    if (use_gw_front_igw) then
      call shr_assert(all(unset_r8 /= [effgw_cm_igw, taubgnd_igw]), &
        'gw_drag_init: Frontogenesis inertial waves enabled, but not all required namelist variables were set!' // &
        errMsg(__FILE__, __LINE__))
      if (masterproc) then
        write(iulog, *) 'gw_init: gw spectrum taubgnd_igw, effgw_cm_igw = ', taubgnd_igw, effgw_cm_igw
        write(iulog, *) ' '
      end if
      cm_igw_desc = gaussian_cm_desc(band_long, kbot_front, kfront, frontgfc, taubgnd_igw, front_gaussian_width)
      ! Output for gravity waves from frontogenesis.
      call gw_spec_addflds(prefix=cm_igw_pf, scheme='C&M IGW', band=band_long, history_defaults=history_waccm)
    end if

    if (use_gw_convect_dp) then
      ttend_dp_idx = pbuf_get_index('TTEND_DP')
      ! Set the deep scheme specification components.
      beres_dp_desc%storm_shift = .true.
      do k = 0, pver
        ! 700 hPa index
        if (pref_edge(k+1) < 70000.0_r8) beres_dp_desc%k = k + 1
      end do
      if (masterproc) then
        write(iulog,*) 'Beres deep level =', beres_dp_desc%k
      end if
      ! Don't use deep convection heating depths below this limit.
      ! This is important for QBO. Bad result if left at 2.5 km.
      beres_dp_desc%min_hdepth = 1000.0_r8
      ! Read Beres file.
      call shr_assert(trim(gw_drag_file) /= '', &
        'gw_drag_init: No gw_drag_file provided for Beres deep scheme. Set this via namelist.' // &
        errMsg(__FILE__, __LINE__))
      call gw_init_beres(gw_drag_file, band_mid, beres_dp_desc)
      ! Output for gravity waves from the Beres scheme (deep).
      call gw_spec_addflds(prefix=beres_dp_pf, scheme='Beres (deep)', &
            band=band_mid, history_defaults=history_waccm)
      call addfld('NETDT' , ['lev'] , 'A', 'K/s'  , 'Net heating rate')
      call addfld('MAXQ0' , horiz_only, 'A', 'K/day', 'Max column heating rate')
      call addfld('HDEPTH', horiz_only, 'A', 'km'   , 'Heating Depth')
      if (history_waccm) then
        call add_default('NETDT    ', 1, ' ')
        call add_default('HDEPTH   ', 1, ' ')
        call add_default('MAXQ0    ', 1, ' ')
      end if
    end if
    if (use_gw_convect_sh) then
      ttend_sh_idx = pbuf_get_index('TTEND_SH')
      ! Set the shallow scheme specification components.
      beres_sh_desc%storm_shift = .false.
      do k = 0, pver
        ! 900 hPa index
        if (pref_edge(k+1) < 90000.0_r8) beres_sh_desc%k = k + 1
      end do
      if (masterproc) then
        write (iulog, *) 'Beres shallow level =', beres_sh_desc%k
      end if
      ! Use all heating depths for shallow convection.
      beres_sh_desc%min_hdepth = 0
      ! Read Beres file.
      call shr_assert(trim(gw_drag_file_sh) /= "", &
            "gw_drag_init: No gw_drag_file_sh provided for Beres shallow &
            &scheme. Set this via namelist."// &
            errMsg(__FILE__, __LINE__))
      call gw_init_beres(gw_drag_file_sh, band_mid, beres_sh_desc)
      ! Output for gravity waves from the Beres scheme (shallow).
      call gw_spec_addflds(prefix=beres_sh_pf, scheme="Beres (shallow)", &
            band=band_mid, history_defaults=history_waccm)
      call addfld('SNETDT' , ['lev'] , 'A', 'K/s'  , 'Net heating rate')
      call addfld('SMAXQ0' , horiz_only, 'A', 'K/day', 'Max column heating rate')
      call addfld('SHDEPTH', horiz_only, 'A', 'km'   , 'Heating Depth')
      if (history_waccm) then
        call add_default('SNETDT   ', 1, ' ')
        call add_default('SHDEPTH  ', 1, ' ')
        call add_default('SMAXQ0   ', 1, ' ')
      end if
    end if

    call addfld('EKGW' ,[ 'ilev' ], 'A','M2/S', &
        'Effective Kzz due to diffusion by gravity waves')

    if (history_waccm) then
      call add_default('EKGW', 1, ' ')
    end if

    call addfld('UTGW_TOTAL',    [ 'lev' ], 'A','m/s2', &
        'Total U tendency due to gravity wave drag')
    call addfld('VTGW_TOTAL',    [ 'lev' ], 'A','m/s2', &
        'Total V tendency due to gravity wave drag')
    call register_vector_field('UTGW_TOTAL', 'VTGW_TOTAL')

    ! Total temperature tendency output.
    call addfld('TTGW', [ 'lev' ], 'A', 'K/s',  &
        'T tendency - gravity wave drag')

    ! Water budget terms.
    call addfld('QTGW',[ 'lev' ], 'A','kg/kg/s', &
        'Q tendency - gravity wave drag')
    call addfld('CLDLIQTGW',[ 'lev' ], 'A','kg/kg/s', &
        'CLDLIQ tendency - gravity wave drag')
    call addfld('CLDICETGW',[ 'lev' ], 'A','kg/kg/s', &
        'CLDICE tendency - gravity wave drag')

    if (history_budget) then
      call add_default('TTGW', history_budget_histfile_num, ' ')
      call add_default('QTGW', history_budget_histfile_num, ' ')
      call add_default('CLDLIQTGW', history_budget_histfile_num, ' ')
      call add_default('CLDICETGW', history_budget_histfile_num, ' ')
    end if

    ! Get indices to actually output the above.
    call cnst_get_ind("CLDLIQ", ixcldliq)
    call cnst_get_ind("CLDICE", ixcldice)

    if (gw_top_taper) then
      allocate(vramp(pver))
      vramp  = 1
      topndx = 1
      botndx = press_lim_idx(0.6e-02_r8, top=.true.)
      if (botndx>1) then
        do k = botndx, topndx, -1
          vramp(k) = vramp(k+1) / (pref_edge(k+1) / pref_edge(k))
        end do
        if (masterproc) then
          write(iulog, '(A)') 'GW taper coef (vramp):'
          do k = 1, pver
            write(iulog, "('k: ',I4,' taper coef,press(Pa): ',F12.8,E12.4)") k, vramp(k), pref_mid(k)
          end do
        end if
      end if
    end if

  end subroutine gw_init

  subroutine gw_init_beres(file_name, band, desc)

    use ioFileMod, only: getfil
    use pio, only: file_desc_t, pio_nowrite, pio_inq_varid, pio_get_var, &
        pio_closefile
    use cam_pio_utils, only: cam_pio_openfile

    character(*), intent(in) :: file_name
    type(GWBand), intent(in) :: band
    type(BeresSourceDesc), intent(inout) :: desc
  
    type(file_desc_t) gw_file_desc
    ! PIO variable ids and error code.
    integer mfccid, hdid, stat
    ! Number of wavenumbers in the input file.
    integer ngwv_file
    ! Full path to gw_drag_file.
    character(256) file_path
    character(256) msg
  
    !----------------------------------------------------------------------
    ! Read in look-up table for source spectra
    !-----------------------------------------------------------------------
  
    call getfil(file_name, file_path)
    call cam_pio_openfile(gw_file_desc, file_path, pio_nowrite)
  
    ! Get HD (heating depth) dimension.
    desc%maxh  = get_pio_dimlen(gw_file_desc, 'HD', file_path)
  
    ! Get MW (mean wind) dimension.
  
    desc%maxuh = get_pio_dimlen(gw_file_desc, 'MW', file_path)
  
    ! Get PS (phase speed) dimension.
  
    ngwv_file  = get_pio_dimlen(gw_file_desc, 'PS', file_path)
  
    ! Number in each direction is half of total (and minus phase speed of 0).
    desc%maxuh = (desc%maxuh - 1) / 2
    ngwv_file  = (ngwv_file - 1) / 2
  
    call shr_assert(ngwv_file >= band%ngwv, &
         'gw_beres_init: PS in lookup table file does not cover the whole &
         &spectrum implied by the model''s ngwv.')
  
    ! Allocate hd and get data.
  
    allocate(desc%hd(desc%maxh), stat=stat, errmsg=msg)
  
    call shr_assert(stat == 0, &
         'gw_init_beres: Allocation error (hd): '//msg// &
         errMsg(__FILE__, __LINE__))
  
    stat = pio_inq_varid(gw_file_desc, 'HD', hdid)
  
    call handle_pio_error(stat, 'Error finding HD in: ' // trim(file_path))
  
    stat = pio_get_var(gw_file_desc, hdid, start=[1], count=[desc%maxh], ival=desc%hd)
  
    call handle_pio_error(stat, 'Error reading HD from: ' // trim(file_path))
  
    ! While not currently documented in the file, it uses kilometers. Convert
    ! to meters.
    desc%hd = desc%hd * 1000.0_r8
  
    ! Allocate mfcc. "desc%maxh" and "desc%maxuh" are from the file, but the
    ! model determines wavenumber dimension.
  
    allocate(desc%mfcc(desc%maxh,-desc%maxuh:desc%maxuh,-band%ngwv:band%ngwv), stat=stat, errmsg=msg)
  
    call shr_assert(stat == 0, &
      'gw_init_beres: Allocation error (mfcc): '//msg// &
      errMsg(__FILE__, __LINE__))
  
    ! Get mfcc data.
    stat = pio_inq_varid(gw_file_desc, 'mfcc', mfccid)
  
    call handle_pio_error(stat, 'Error finding mfcc in: ' // trim(file_path))
  
    stat = pio_get_var(gw_file_desc, mfccid, &
      start=[1,1,ngwv_file-band%ngwv+1], count=shape(desc%mfcc), &
      ival=desc%mfcc)
  
    call handle_pio_error(stat, 'Error reading mfcc from: ' // trim(file_path))
  
    call pio_closefile(gw_file_desc)
  
    if (masterproc) then
      write(iulog, *) 'Read in source spectra from file.'
      write(iulog, *) 'mfcc max, min = ', maxval(desc%mfcc), ', ', minval(desc%mfcc)
    end if
  
  end subroutine gw_init_beres

  ! Utility to reduce the repetitiveness of reads during initialization.
  function get_pio_dimlen(file_desc, dim_name, file_path) result(dimlen)

    use pio, only: file_desc_t, pio_inq_dimid, pio_inq_dimlen

    type(file_desc_t), intent(in) :: file_desc
    character(*), intent(in) :: dim_name

    ! File path, for use in error messages only.
    character(*), intent(in) :: file_path

    integer :: dimlen
    integer :: dimid, stat

    stat = pio_inq_dimid(file_desc, dim_name, dimid)

    call handle_pio_error(stat, &
       'Error finding dimension '//dim_name//' in: '//file_path)

    stat = pio_inq_dimlen(file_desc, dimid, dimlen)

    call handle_pio_error(stat, &
       'Error reading dimension '//dim_name//' from: '//file_path)

  end function get_pio_dimlen

  ! In fact, we'd usually expect PIO errors to abort the run before you can
  ! even check the error code. But just in case, use this little assert.
  subroutine handle_pio_error(stat, message)

    use pio, only: pio_noerr

    integer, intent(in) :: stat
    character(*) :: message

    call shr_assert(stat == pio_noerr, &
       'PIO error in gw_init_beres: '//trim(message)// &
       errMsg(__FILE__, __LINE__))

  end subroutine handle_pio_error

  subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)

    use physics_types,  only: physics_state_copy, set_dry_to_wet
    use constituents,   only: cnst_type
    use physics_buffer, only: physics_buffer_desc, pbuf_get_field
    use camsrfexch, only: cam_in_t
    ! Location-dependent cpair
    use physconst,  only: cpairv, pi
    use coords_1d,  only: Coords1D
    use gw_common,  only: gw_prof, gw_drag_prof, calc_taucd, &
         momentum_flux, momentum_fixer, energy_change, energy_fixer, &
         coriolis_speed, adjust_inertial
    use gw_oro,     only: gw_oro_src
    use gw_front,   only: gw_cm_src
    use gw_convect, only: gw_beres_src

    type(physics_state), intent(in) :: state   ! physics state structure
    type(physics_buffer_desc), pointer :: pbuf(:) ! Physics buffer
    real(r8), intent(in) :: dt                    ! time step
    ! Parameterization net tendencies.
    type(physics_ptend), intent(out):: ptend
    type(cam_in_t), intent(in) :: cam_in
    real(r8), intent(out) :: flx_heat(pcols)

    type(physics_state) :: state1     ! Local copy of state variable

    integer lchnk                  ! chunk identifier
    integer ncol                   ! number of atmospheric columns
    integer i, k                   ! loop indices
    type(Coords1D) :: p               ! Pressure coordinates

    real(r8) ttgw(state%ncol,pver) ! temperature tendency
    real(r8) utgw(state%ncol,pver) ! zonal wind tendency
    real(r8) vtgw(state%ncol,pver) ! meridional wind tendency

    real(r8) ni(state%ncol,pver+1) ! interface Brunt-Vaisalla frequency
    real(r8) nm(state%ncol,pver)   ! midpoint Brunt-Vaisalla frequency
    real(r8) rhoi(state%ncol,pver+1)     ! interface density
    real(r8), allocatable :: tau(:,:,:)  ! wave Reynolds stress
    real(r8) tau0x(state%ncol)     ! c=0 sfc. stress (zonal)
    real(r8) tau0y(state%ncol)     ! c=0 sfc. stress (meridional)
    real(r8) ubi(state%ncol,pver+1)! projection of wind at interfaces
    real(r8) ubm(state%ncol,pver)  ! projection of wind at midpoints
    real(r8) xv(state%ncol)        ! unit vector of source wind (x)
    real(r8) yv(state%ncol)        ! unit vector of source wind (y)

    integer m
    real(r8) qtgw(state%ncol,pver,pcnst) ! constituents tendencies

    ! Reynolds stress for waves propagating in each cardinal direction.
    real(r8) taucd(state%ncol,pver+1,4)

    ! gravity wave wind tendency for each wave
    real(r8), allocatable :: gwut(:,:,:)

    ! Temperature tendencies from diffusion and kinetic energy.
    real(r8) dttdf(state%ncol,pver)
    real(r8) dttke(state%ncol,pver)

    ! Wave phase speeds for each column
    real(r8), allocatable :: c(:,:)

    ! Efficiency for a gravity wave source.
    real(r8) effgw(state%ncol)

    ! Coriolis characteristic speed.
    real(r8) u_coriolis(state%ncol)

    ! Adjustment for inertial gravity waves.
    real(r8), allocatable :: ro_adjust(:,:,:)

    ! pbuf fields
    ! Molecular diffusivity
    real(r8), pointer :: kvt_in(:,:)
    real(r8) kvtt(state%ncol,pver+1)

    ! Frontogenesis
    real(r8), pointer :: frontgf(:,:)
    real(r8), pointer :: frontga(:,:)

    ! Temperature change due to deep convection.
    real(r8), pointer :: ttend_dp(:,:)
    ! Temperature change due to shallow convection.
    real(r8), pointer :: ttend_sh(:,:)

    ! Standard deviation of orography.
    real(r8), pointer :: sgh(:)

    ! gridbox area
    real(r8), pointer :: gbxar(:)

       ! Beta ridges
    ! width of ridges.
    real(r8), pointer :: hwdth(:,:)
    ! length of ridges.
    real(r8), pointer :: clngt(:,:)
    ! Maximum deviations of ridges.
    real(r8), pointer :: mxdis(:,:)
    ! orientation of ridges.
    real(r8), pointer :: angll(:,:)
    ! anisotropy of ridges.
    real(r8), pointer :: anixy(:,:)

       ! Gamma ridges
    ! width of ridges.
    real(r8), pointer :: hwdthg(:,:)
    ! length of ridges.
    real(r8), pointer :: clngtg(:,:)
    ! Maximum deviations of ridges.
    real(r8), pointer :: mxdisg(:,:)
    ! orientation of ridges.
    real(r8), pointer :: angllg(:,:)
    ! anisotropy of ridges.
    real(r8), pointer :: anixyg(:,:)

    ! Indices of gravity wave source and lowest level where wind tendencies
    ! are allowed.
    integer src_level(state%ncol)
    integer tend_level(state%ncol)

    ! Convective source heating depth.
    ! heating depth
    real(r8) hdepth(state%ncol)
    ! maximum heating rate
    real(r8) maxq0(state%ncol)

    ! Scale sgh to account for landfrac.
    real(r8) sgh_scaled(state%ncol)

    ! Parameters for the IGW polar taper.
    real(r8), parameter :: degree2radian = pi/180.0_r8
    real(r8), parameter :: al0 = 82.5_r8 * degree2radian
    real(r8), parameter :: dlat0 = 5.0_r8 * degree2radian

    ! effective gw diffusivity at interfaces needed for output
    real(r8) egwdffi(state%ncol,pver+1)
    ! sum from the two types of spectral GW
    real(r8) egwdffi_tot(state%ncol,pver+1)

    ! Momentum fluxes used by fixer.
    real(r8) um_flux(state%ncol), vm_flux(state%ncol)
    ! Energy change used by fixer.
    real(r8) de(state%ncol)

    ! Which constituents are being affected by diffusion.
    logical  lq(pcnst)

    ! Contiguous copies of state arrays.
    real(r8) dse (state%ncol,pver)
    real(r8) t   (state%ncol,pver)
    real(r8) u   (state%ncol,pver)
    real(r8) v   (state%ncol,pver)
    real(r8) q   (state%ncol,pver,pcnst)
    real(r8) piln(state%ncol,pver+1)
    real(r8) zm  (state%ncol,pver)
    real(r8) zi  (state%ncol,pver+1)

    ! Make local copy of input state.
    call physics_state_copy(state, state1)

    ! constituents are all treated as wet mmr
    call set_dry_to_wet(state1)

    lchnk = state1%lchnk
    ncol  = state1%ncol
    p     = Coords1D(state1%pint(:ncol,:))
    dse   = state1%s     (:ncol,:)
    t     = state1%t     (:ncol,:)
    u     = state1%u     (:ncol,:)
    v     = state1%v     (:ncol,:)
    q     = state1%q     (:ncol,:,:)
    piln  = state1%lnpint(:ncol,:)
    zm    = state1%zm    (:ncol,:)
    zi    = state1%zi    (:ncol,:)

    lq = .true.
    call physics_ptend_init(ptend, state1%psetcols, 'Gravity wave drag', &
         ls=.true., lu=.true., lv=.true., lq=lq)

    ! Profiles of background state variables
    call gw_prof(ncol, p, cpair, t, rhoi, nm, ni)

    if (do_molec_diff) then
      !--------------------------------------------------------
      ! Initialize and calculate local molecular diffusivity
      !--------------------------------------------------------

      call pbuf_get_field(pbuf, kvt_idx, kvt_in)  ! kvt_in(1:pcols,1:pver+1)

      ! Set kvtt from pbuf field; kvtt still needs a factor of 1/cpairv.
      kvtt = kvt_in(:ncol,:)

      ! Use linear extrapolation of cpairv to top interface.
      kvtt(:,1) = kvtt(:,1) / &
           (1.5_r8*cpairv(:ncol,1,lchnk) - &
           0.5_r8*cpairv(:ncol,2,lchnk))

      ! Interpolate cpairv to other interfaces.
      do k = 2, nbot_molec
        kvtt(:,k) = kvtt(:,k) / &
              (cpairv(:ncol,k+1,lchnk)+cpairv(:ncol,k,lchnk)) * 2._r8
      end do
    else
      kvtt = 0.0_r8
    end if

    if (use_gw_front_igw) then
      u_coriolis = coriolis_speed(band_long, state1%lat(:ncol))
    end if

    ! Totals that accumulate over different sources.
    egwdffi_tot = 0.0_r8
    flx_heat = 0.0_r8

    if (use_gw_convect_dp) then
      !------------------------------------------------------------------
      ! Convective gravity waves (Beres scheme, deep).
      !------------------------------------------------------------------

      ! Allocate wavenumber fields.
      allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1))
      allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv))
      allocate(c(ncol,-band_mid%ngwv:band_mid%ngwv))

      ! Set up heating
      call pbuf_get_field(pbuf, ttend_dp_idx, ttend_dp)

      ! Efficiency of gravity wave momentum transfer.
      ! This is really only to remove the pole points.
      where (pi/2._r8 - abs(state1%lat(:ncol)) >= 4*epsilon(1._r8))
        effgw = effgw_beres_dp
      else where
        effgw = 0.0_r8
      end where

      ! Determine wave sources for Beres deep scheme
      call gw_beres_src(ncol, band_mid, beres_dp_desc, &
           u, v, ttend_dp(:ncol,:), zm, src_level, tend_level, tau, &
           ubm, ubi, xv, yv, c, hdepth, maxq0)

      ! Solve for the drag profile with Beres source spectrum.
      call gw_drag_prof(ncol, band_mid, p, src_level, tend_level, dt, &
           t, vramp,    &
           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
           effgw,   c,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
           lapply_effgw_in=gw_apply_tndmax)

      ! Project stress into directional components.
      taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, c, xv, yv, ubi)

      !  add the diffusion coefficients
      do k = 1, pver+1
        egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
      end do

      ! Store constituents tendencies
      do m = 1, pcnst
        do k = 1, pver
          ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        end do
      end do

      ! Find momentum flux, and use it to fix the wind tendencies below
      ! the gravity wave region.
      call momentum_flux(tend_level, taucd, um_flux, vm_flux)
      call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)

      ! Add the momentum tendencies to the output tendency arrays.
      do k = 1, pver
        ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
      end do

      ! Find energy change in the current state, and use fixer to apply
      ! the difference in lower levels.
      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
           ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
      call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)

      do k = 1, pver
        ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
      end do

      ! Change ttgw to a temperature tendency before outputing it.
      ttgw = ttgw / cpair
      call gw_spec_outflds(beres_dp_pf, lchnk, ncol, band_mid, c, u, v, &
           xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
           taucd)

      ! Diagnostic outputs (convert hdepth to km).
      call outfld('NETDT', ttend_dp, pcols, lchnk)
      call outfld('HDEPTH', hdepth/1000.0_r8, ncol, lchnk)
      call outfld('MAXQ0', maxq0, ncol, lchnk)

      deallocate(tau, gwut, c)
    end if

    if (use_gw_convect_sh) then
      !------------------------------------------------------------------
      ! Convective gravity waves (Beres scheme, shallow).
      !------------------------------------------------------------------

      ! Allocate wavenumber fields.
      allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1))
      allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv))
      allocate(c(ncol,-band_mid%ngwv:band_mid%ngwv))

      ! Set up heating
      call pbuf_get_field(pbuf, ttend_sh_idx, ttend_sh)

      ! Efficiency of gravity wave momentum transfer.
      ! This is really only to remove the pole points.
      where (pi/2._r8 - abs(state1%lat(:ncol)) >= 4*epsilon(1._r8))
        effgw = effgw_beres_sh
      else where
        effgw = 0.0_r8
      end where

      ! Determine wave sources for Beres shallow scheme
      call gw_beres_src(ncol, band_mid, beres_sh_desc, &
           u, v, ttend_sh(:ncol,:), zm, src_level, tend_level, tau, &
           ubm, ubi, xv, yv, c, hdepth, maxq0)

      ! Solve for the drag profile with Beres source spectrum.
      call gw_drag_prof(ncol, band_mid, p, src_level, tend_level,  dt, &
           t, vramp,    &
           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
           effgw,   c,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
           lapply_effgw_in=gw_apply_tndmax)

      ! Project stress into directional components.
      taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, c, xv, yv, ubi)

      !  add the diffusion coefficients
      do k = 1, pver + 1
        egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
      end do

      ! Store constituents tendencies
      do m = 1, pcnst
        do k = 1, pver
          ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        end do
      end do

      ! Add the momentum tendencies to the output tendency arrays.
      ! Don't calculate fixers, since we are too close to the ground to
      ! spread momentum/energy differences across low layers.
      do k = 1, pver
        ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
      end do

      ! Calculate energy change for output to CAM's energy checker.
      ! This is sort of cheating; we don't have a good a priori idea of the
      ! energy coming from surface stress, so we just integrate what we and
      ! actually have so far and overwrite flx_heat with that.
      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
           ptend%v(:ncol,:), ptend%s(:ncol,:), de)
      flx_heat(:ncol) = de

      ! Change ttgw to a temperature tendency before outputing it.
      ttgw = ttgw / cpair
      call gw_spec_outflds(beres_sh_pf, lchnk, ncol, band_mid, c, u, v, &
           xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
           taucd)

      ! Diagnostic outputs (convert SHDEPTH to km).
      call outfld ('SNETDT', ttend_sh, pcols, lchnk)
      call outfld ('SHDEPTH', hdepth/1000.0_r8, ncol, lchnk)
      call outfld ('SMAXQ0', maxq0, ncol, lchnk)

      deallocate(tau, gwut, c)
    end if

    if (use_gw_front .or. use_gw_front_igw) then
      ! Get frontogenesis physics buffer fields set by dynamics.
      call pbuf_get_field(pbuf, frontgf_idx, frontgf)
      call pbuf_get_field(pbuf, frontga_idx, frontga)

      ! Output for diagnostics.
      call outfld ('FRONTGF', frontgf, pcols, lchnk)
      call outfld ('FRONTGFA', frontga, pcols, lchnk)
    end if

    if (use_gw_front) then
      !------------------------------------------------------------------
      ! Frontally generated gravity waves
      !------------------------------------------------------------------

      ! Allocate wavenumber fields.
      allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1))
      allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv))
      allocate(c(ncol,-band_mid%ngwv:band_mid%ngwv))

      ! Efficiency of gravity wave momentum transfer.
      effgw = effgw_cm
      ! Frontogenesis is too high at the poles (at least for the FV
      ! dycore), so introduce a polar taper.
      if (gw_polar_taper) effgw = effgw * cos(state1%lat(:ncol))

      ! Determine the wave source for C&M background spectrum
      call gw_cm_src(ncol, band_mid, cm_desc, u, v, frontgf(:ncol,:), &
           src_level, tend_level, tau, ubm, ubi, xv, yv, c)

      ! Solve for the drag profile with C&M source spectrum.
      call gw_drag_prof(ncol, band_mid, p, src_level, tend_level,  dt, &
           t, vramp,   &
           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
           effgw,   c,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
           lapply_effgw_in=gw_apply_tndmax)

      ! Project stress into directional components.
      taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, c, xv, yv, ubi)

      !  add the diffusion coefficients
      do k = 1, pver+1
        egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
      end do

      !Add the constituent tendencies
      do m=1, pcnst
        do k = 1, pver
          ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        end do
      end do

      ! Find momentum flux, and use it to fix the wind tendencies below
      ! the gravity wave region.
      call momentum_flux(tend_level, taucd, um_flux, vm_flux)
      call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)

      ! add the momentum tendencies to the output tendency arrays
      do k = 1, pver
        ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
      end do

      ! Find energy change in the current state, and use fixer to apply
      ! the difference in lower levels.
      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
            ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
      call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)

      do k = 1, pver
        ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
      end do

      ! Change ttgw to a temperature tendency before outputing it.
      ttgw = ttgw / cpair
      call gw_spec_outflds(cm_pf, lchnk, ncol, band_mid, c, u, v, &
            xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
            taucd)
      deallocate(tau, gwut, c)
    end if

    if (use_gw_front_igw) then
      !------------------------------------------------------------------
      ! Frontally generated inertial gravity waves
      !------------------------------------------------------------------

      ! Allocate wavenumber fields.
      allocate(tau(ncol,-band_long%ngwv:band_long%ngwv,pver+1))
      allocate(gwut(ncol,pver,-band_long%ngwv:band_long%ngwv))
      allocate(c(ncol,-band_long%ngwv:band_long%ngwv))
      allocate(ro_adjust(ncol,-band_long%ngwv:band_long%ngwv,pver+1))

      ! Efficiency of gravity wave momentum transfer.
      effgw = effgw_cm_igw

      ! Frontogenesis is too high at the poles (at least for the FV
      ! dycore), so introduce a polar taper.
      if (gw_polar_taper) then
        where (abs(state1%lat(:ncol)) <= 89._r8*degree2radian)
          effgw = effgw * 0.25_r8 * &
              (1._r8+tanh((state1%lat(:ncol)+al0)/dlat0)) * &
              (1._r8-tanh((state1%lat(:ncol)-al0)/dlat0))
        else where
          effgw = 0.0_r8
        end where
      end if

      ! Determine the wave source for C&M background spectrum
      call gw_cm_src(ncol, band_long, cm_igw_desc, u, v, frontgf(:ncol,:), &
            src_level, tend_level, tau, ubm, ubi, xv, yv, c)

      call adjust_inertial(band_long, tend_level, u_coriolis, c, ubi, &
            tau, ro_adjust)

      ! Solve for the drag profile with C&M source spectrum.
      call gw_drag_prof(ncol, band_long, p, src_level, tend_level,  dt, &
            t, vramp,    &
            piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
            effgw,   c,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
            ttgw, qtgw, egwdffi,  gwut, dttdf, dttke, ro_adjust=ro_adjust, &
            lapply_effgw_in=gw_apply_tndmax)

      ! Project stress into directional components.
      taucd = calc_taucd(ncol, band_long%ngwv, tend_level, tau, c, xv, yv, ubi)

      ! Add the diffusion coefficients
      do k = 1, pver+1
        egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
      end do

      ! Add the constituent tendencies
      do m=1, pcnst
        do k = 1, pver
          ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        end do
      end do

      ! Find momentum flux, and use it to fix the wind tendencies below
      ! the gravity wave region.
      call momentum_flux(tend_level, taucd, um_flux, vm_flux)
      call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)

      ! add the momentum tendencies to the output tendency arrays
      do k = 1, pver
        ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
      end do

      ! Find energy change in the current state, and use fixer to apply
      ! the difference in lower levels.
      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
            ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
      call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)

      do k = 1, pver
        ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
      end do

      ! Change ttgw to a temperature tendency before outputing it.
      ttgw = ttgw / cpair
      call gw_spec_outflds(cm_igw_pf, lchnk, ncol, band_long, c, u, v, &
            xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
            taucd)
      deallocate(tau, gwut, c, ro_adjust)
    end if

    if (use_gw_oro) then
      !---------------------------------------------------------------------
      ! Orographic stationary gravity waves
      !---------------------------------------------------------------------

      ! Allocate wavenumber fields.
      allocate(tau(ncol,band_oro%ngwv:band_oro%ngwv,pver+1))
      allocate(gwut(ncol,pver,band_oro%ngwv:band_oro%ngwv))
      allocate(c(ncol,band_oro%ngwv:band_oro%ngwv))

      ! Efficiency of gravity wave momentum transfer.
      ! Take into account that wave sources are only over land.
      call pbuf_get_field(pbuf, sgh_idx, sgh)

      if (gw_lndscl_sgh) then
        where (cam_in%landfrac(:ncol) >= epsilon(1._r8))
          effgw = effgw_oro * cam_in%landfrac(:ncol)
          sgh_scaled = sgh(:ncol) / sqrt(cam_in%landfrac(:ncol))
        else where
          effgw = 0.0_r8
          sgh_scaled = 0.0_r8
        end where

        ! Determine the orographic wave source
        call gw_oro_src(ncol, band_oro, p, &
          u, v, t, sgh_scaled, zm, nm, &
          src_level, tend_level, tau, ubm, ubi, xv, yv, c)
      else
        effgw = effgw_oro

        ! Determine the orographic wave source
        call gw_oro_src(ncol, band_oro, p, &
          u, v, t, sgh(:ncol), zm, nm, &
          src_level, tend_level, tau, ubm, ubi, xv, yv, c)
      endif
      do i = 1, ncol
        if (state1%lat(i) < 0.0_r8) then
          tau(i,:,:) = tau(i,:,:) * gw_oro_south_fac
        end if
      end do

      ! Solve for the drag profile with orographic sources.
      call gw_drag_prof(ncol, band_oro, p, src_level, tend_level,   dt,   &
           t, vramp,   &
           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
           effgw,c,          kvtt, q,  dse,  tau,  utgw,  vtgw, &
           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
           lapply_effgw_in=gw_apply_tndmax)

      ! For orographic waves, don't bother with taucd, since there are no
      ! momentum conservation routines or directional diagnostics.

      !  add the diffusion coefficients
      do k = 1, pver+1
        egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
      end do

      ! Add the orographic tendencies to the spectrum tendencies.
      ! Don't calculate fixers, since we are too close to the ground to
      ! spread momentum/energy differences across low layers.
      do k = 1, pver
        ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
        ! Convert to temperature tendency for output.
        ttgw(:,k) = ttgw(:,k) / cpairv(:ncol, k, lchnk)
      end do

      ! Calculate energy change for output to CAM's energy checker.
      ! This is sort of cheating; we don't have a good a priori idea of the
      ! energy coming from surface stress, so we just integrate what we and
      ! actually have so far and overwrite flx_heat with that.
      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
           ptend%v(:ncol,:), ptend%s(:ncol,:), de)
      flx_heat(:ncol) = de

      do m = 1, pcnst
        do k = 1, pver
          ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        end do
      end do

      ! Write output fields to history file
      call outfld('TAUAORO', tau(:,0,:),  ncol, lchnk)
      call outfld('UTGWORO', utgw,  ncol, lchnk)
      call outfld('VTGWORO', vtgw,  ncol, lchnk)
      call outfld('TTGWORO', ttgw,  ncol, lchnk)
      call outfld('TTGWSDFORO', dttdf / cpair,  ncol, lchnk)
      call outfld('TTGWSKEORO', dttke / cpair,  ncol, lchnk)
      tau0x = tau(:,0,pver+1) * xv
      tau0y = tau(:,0,pver+1) * yv
      call outfld('TAUGWX', tau0x, ncol, lchnk)
      call outfld('TAUGWY', tau0y, ncol, lchnk)

      deallocate(tau, gwut, c)
    end if

    if (use_gw_rdg_beta) then
      !---------------------------------------------------------------------
      ! Orographic stationary gravity waves
      !---------------------------------------------------------------------

      ! Efficiency of gravity wave momentum transfer.
      ! Take into account that wave sources are only over land.

      hwdth => rdg_hwdth(:ncol,:,lchnk)
      clngt => rdg_clngt(:ncol,:,lchnk)
      gbxar => rdg_gbxar(:ncol,lchnk)
      mxdis => rdg_mxdis(:ncol,:,lchnk)
      angll => rdg_angll(:ncol,:,lchnk)
      anixy => rdg_anixy(:ncol,:,lchnk)

      where (mxdis < 0.0_r8)
        mxdis = 0.0_r8
      end where

      ! Save state at top of routine
      ! Useful for unit testing checks
      call outfld('UEGW', u ,  ncol, lchnk)
      call outfld('VEGW', v ,  ncol, lchnk)
      call outfld('TEGW', t ,  ncol, lchnk)

      call gw_rdg_calc(&
         'BETA ', ncol, lchnk, n_rdg_beta, dt,     &
         u, v, t, p, piln, zm, zi,                 &
         nm, ni, rhoi, kvtt, q, dse,               &
         effgw_rdg_beta, effgw_rdg_beta_max,       &
         hwdth, clngt, gbxar, mxdis, angll, anixy, &
         rdg_beta_cd_llb, trpd_leewv_rdg_beta,     &
         ptend, flx_heat)
    end if

    if (use_gw_rdg_gamma) then
      !---------------------------------------------------------------------
      ! Orographic stationary gravity waves
      !---------------------------------------------------------------------

      ! Efficiency of gravity wave momentum transfer.
      ! Take into account that wave sources are only over land.

      hwdthg => rdg_hwdthg(:ncol,:,lchnk)
      clngtg => rdg_clngtg(:ncol,:,lchnk)
      gbxar  => rdg_gbxar(:ncol,lchnk)
      mxdisg => rdg_mxdisg(:ncol,:,lchnk)
      angllg => rdg_angllg(:ncol,:,lchnk)
      anixyg => rdg_anixyg(:ncol,:,lchnk)

      where(mxdisg < 0.0_r8)
        mxdisg = 0.0_r8
      end where

      call gw_rdg_calc(&
          'GAMMA', ncol, lchnk, n_rdg_gamma, dt,         &
          u, v, t, p, piln, zm, zi,                      &
          nm, ni, rhoi, kvtt, q, dse,                    &
          effgw_rdg_gamma, effgw_rdg_gamma_max,          &
          hwdthg, clngtg, gbxar, mxdisg, angllg, anixyg, &
          rdg_gamma_cd_llb, trpd_leewv_rdg_gamma,        &
          ptend, flx_heat)
    end if

    ! Convert the tendencies for the dry constituents to dry air basis.
    do m = 1, pcnst
      if (cnst_type(m) == 'dry') then
        do k = 1, pver
          do i = 1, ncol
            ptend%q(i,k,m) = ptend%q(i,k,m)*state1%pdel(i,k)/state1%pdeldry(i,k)
          end do
        end do
      end if
    end do

    ! Write totals to history file.
    call outfld('EKGW', egwdffi_tot , ncol, lchnk)
    call outfld('TTGW', ptend%s/cpairv(:,:,lchnk),  pcols, lchnk)

    call outfld('UTGW_TOTAL', ptend%u, pcols, lchnk)
    call outfld('VTGW_TOTAL', ptend%v, pcols, lchnk)

    call outfld('QTGW', ptend%q(:,:,1), pcols, lchnk)
    call outfld('CLDLIQTGW', ptend%q(:,:,ixcldliq), pcols, lchnk)
    call outfld('CLDICETGW', ptend%q(:,:,ixcldice), pcols, lchnk)

    ! Destroy objects.
    call p%finalize()

  end subroutine gw_tend

  subroutine gw_rdg_calc( &
    type, ncol, lchnk, n_rdg, dt, &
    u, v, t, p, piln, zm, zi, &
    nm, ni, rhoi, kvtt, q, dse, &
    effgw_rdg, effgw_rdg_max, &
    hwdth, clngt, gbxar, &
    mxdis, angll, anixy, &
    rdg_cd_llb, trpd_leewv, &
    ptend, flx_heat)

    use coords_1d,  only: Coords1D
    use gw_rdg,     only: gw_rdg_src, gw_rdg_belowpeak, gw_rdg_break_trap, gw_rdg_do_vdiff
    use gw_common,  only: gw_drag_prof, energy_change

    character(5),   intent(in) :: type                ! BETA or GAMMA
    integer ,       intent(in) :: ncol
    integer ,       intent(in) :: lchnk
    integer ,       intent(in) :: n_rdg
    real(r8),       intent(in) :: dt

    real(r8),       intent(in) :: u(ncol,pver)
    real(r8),       intent(in) :: v(ncol,pver)
    real(r8),       intent(in) :: t(ncol,pver)
    type(Coords1D), intent(in) :: p
    real(r8),       intent(in) :: piln(ncol,pver+1)   ! Log of interface pressures.
    real(r8),       intent(in) :: zm(ncol,pver)
    real(r8),       intent(in) :: zi(ncol,pver+1)
    real(r8),       intent(in) :: nm(ncol,pver)       ! Midpoint Brunt-Vaisalla frequencies (s-1).
    real(r8),       intent(in) :: ni(ncol,pver+1)     ! Interface Brunt-Vaisalla frequencies (s-1).
    real(r8),       intent(in) :: rhoi(ncol,pver+1)   ! Interface density (kg m-3).
    real(r8),       intent(in) :: kvtt(ncol,pver+1)   ! Molecular thermal diffusivity.
    real(r8),       intent(in) :: q(:,:,:)
    real(r8),       intent(in) :: dse(ncol,pver)      ! Dry static energy.


    real(r8),       intent(in) :: effgw_rdg           ! Tendency efficiency.
    real(r8),       intent(in) :: effgw_rdg_max
    real(r8),       intent(in) :: hwdth(ncol,prdg)    ! width of ridges.
    real(r8),       intent(in) :: clngt(ncol,prdg)    ! length of ridges.
    real(r8),       intent(in) :: gbxar(ncol)         ! gridbox area

    real(r8),       intent(in) :: mxdis(ncol,prdg)    ! Height estimate for ridge (m).
    real(r8),       intent(in) :: angll(ncol,prdg)    ! orientation of ridges.
    real(r8),       intent(in) :: anixy(ncol,prdg)    ! Anisotropy parameter.

    real(r8),       intent(in) :: rdg_cd_llb          ! Drag coefficient for low-level flow
    logical,        intent(in) :: trpd_leewv

    type(physics_ptend), intent(inout):: ptend

    real(r8),       intent(out) :: flx_heat(pcols)

    integer k, m, nn

    ! Wave Reynolds stress
    real(r8), allocatable :: tau(:,:,:)
    ! Gravity wave wind tendency for each wave
    real(r8), allocatable :: gwut(:,:,:)
    ! Wave phase speeds for each column
    real(r8), allocatable :: c(:,:)

    ! Isotropic source flag [anisotropic orography].
    integer  isoflag(ncol)
    ! Horizontal wavenumber [anisotropic orography].
    real(r8) kwvrdg(ncol)
    ! Efficiency for a gravity wave source.
    real(r8) effgw(ncol)

    ! Indices of top gravity wave source level and lowest level where wind
    ! tendencies are allowed.
    integer  src_level(ncol)
    integer  tend_level(ncol)
    integer  bwv_level(ncol)
    integer  tlb_level(ncol)

    ! Projection of wind at midpoints and interfaces.
    real(r8) ubm(ncol,pver)
    real(r8) ubi(ncol,pver+1)

    ! Unit vectors of source wind (zonal and meridional components).
    real(r8) xv(ncol)
    real(r8) yv(ncol)

    ! Averages over source region.
    real(r8) ubmsrc(ncol) ! On-ridge wind.
    real(r8) usrc(ncol)   ! Zonal wind.
    real(r8) vsrc(ncol)   ! Meridional wind.
    real(r8) nsrc(ncol)   ! B-V frequency.
    real(r8) rsrc(ncol)   ! Density.

    ! normalized wavenumber
    real(r8) m2src(ncol)

    ! Top of low-level flow layer.
    real(r8) tlb(ncol)

    ! Bottom of linear wave region.
    real(r8) bwv(ncol)

    ! Froude numbers for flow/drag regimes
    real(r8) Fr1(ncol)
    real(r8) Fr2(ncol)
    real(r8) Frx(ncol)

    ! Wave Reynolds stresses at source level
    real(r8) tauoro(ncol)
    real(r8) taudsw(ncol)

    ! Surface streamline displacement height for linear waves.
    real(r8) hdspwv(ncol)

    ! Surface streamline displacement height for downslope wind regime.
    real(r8) hdspdw(ncol)

    ! Wave breaking level
    real(r8) wbr(ncol)

    real(r8) utgw(ncol,pver)       ! zonal wind tendency
    real(r8) vtgw(ncol,pver)       ! meridional wind tendency
    real(r8) ttgw(ncol,pver)       ! temperature tendency
    real(r8) qtgw(ncol,pver,pcnst) ! constituents tendencies

    ! Effective gravity wave diffusivity at interfaces.
    real(r8) egwdffi(ncol,pver+1)

    ! Temperature tendencies from diffusion and kinetic energy.
    real(r8) dttdf(ncol,pver)
    real(r8) dttke(ncol,pver)

    ! Wave stress in zonal/meridional direction
    real(r8) taurx(ncol,pver+1)
    real(r8) taurx0(ncol,pver+1)
    real(r8) taury(ncol,pver+1)
    real(r8) taury0(ncol,pver+1)

    ! U,V tendency accumulators
    real(r8) utrdg(ncol,pver)
    real(r8) vtrdg(ncol,pver)
    real(r8) ttrdg(ncol,pver)

    ! Energy change used by fixer.
    real(r8) de(ncol)

    character(1) cn
    character(9) fname(4)

    allocate(tau (ncol,     band_oro%ngwv:band_oro%ngwv,pver+1))
    allocate(gwut(ncol,pver,band_oro%ngwv:band_oro%ngwv))
    allocate(c   (ncol,     band_oro%ngwv:band_oro%ngwv))

    ! Initialize accumulated momentum fluxes and tendencies
    taurx = 0
    taury = 0
    ttrdg = 0
    utrdg = 0
    vtrdg = 0

    do nn = 1, n_rdg
      kwvrdg  = 0.001_r8 / ( hwdth(:,nn) + 0.001_r8 ) ! this cant be done every time step !!!
      isoflag = 0
      effgw   = effgw_rdg * ( hwdth(1:ncol,nn)* clngt(1:ncol,nn) ) / gbxar(1:ncol)
      effgw   = min( effgw_rdg_max , effgw )

      call gw_rdg_src(ncol, band_oro, p, &
         u, v, t, mxdis(:,nn), angll(:,nn), anixy(:,nn), kwvrdg, isoflag, zi, nm, &
         src_level, tend_level, bwv_level, tlb_level, tau, ubm, ubi, xv, yv,  &
         ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, c)

      call gw_rdg_belowpeak(ncol, band_oro, rdg_cd_llb, &
         t, mxdis(:,nn), anixy(:,nn), kwvrdg, &
         zi, nm, ni, rhoi, &
         src_level, tau, &
         ubmsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, &
         tauoro, taudsw, hdspwv, hdspdw)

      call gw_rdg_break_trap(ncol, band_oro, &
         zi, nm, ni, ubm, ubi, rhoi, kwvrdg , bwv, tlb, wbr, &
         src_level, tlb_level, hdspwv, hdspdw,  mxdis(:,nn), &
         tauoro, taudsw, tau, &
         ldo_trapped_waves=trpd_leewv)

      call gw_drag_prof(ncol, band_oro, p, src_level, tend_level, dt, &
         t, vramp,    &
         piln, rhoi, nm, ni, ubm, ubi, xv, yv,   &
         effgw, c, kvtt, q, dse, tau, utgw, vtgw, &
         ttgw, qtgw, egwdffi,   gwut, dttdf, dttke, &
         kwvrdg=kwvrdg, &
         satfac_in = 1._r8, lapply_vdiff=gw_rdg_do_vdiff )

      ! Add the tendencies from each ridge to the totals.
      do k = 1, pver
        ! diagnostics
        utrdg(:,k) = utrdg(:,k) + utgw(:,k)
        vtrdg(:,k) = vtrdg(:,k) + vtgw(:,k)
        ttrdg(:,k) = ttrdg(:,k) + ttgw(:,k)
        ! physics tendencies
        ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
      end do

      do m = 1, pcnst
        do k = 1, pver
          ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        end do
      end do

      do k = 1, pver+1
        taurx0(:,k) =  tau(:,0,k)*xv
        taury0(:,k) =  tau(:,0,k)*yv
        taurx(:,k)  =  taurx(:,k) + taurx0(:,k)
        taury(:,k)  =  taury(:,k) + taury0(:,k)
      end do

      if (nn == 1) then
        call outfld('TAU1RDG'//trim(type)//'M', tau(:,0,:),  ncol, lchnk)
        call outfld('UBM1'//trim(type),         ubm,         ncol, lchnk)
        call outfld('UBT1RDG'//trim(type),      gwut,        ncol, lchnk)
      end if

      if (nn <= 6) then
        write(cn, '(i1)') nn
        call outfld('TAU'//cn//'RDG'//trim(type)//'X', taurx0,  ncol, lchnk)
        call outfld('TAU'//cn//'RDG'//trim(type)//'Y', taury0,  ncol, lchnk)
        call outfld('UT'//cn//'RDG'//trim(type),       utgw,    ncol, lchnk)
        call outfld('VT'//cn//'RDG'//trim(type),       vtgw,    ncol, lchnk)
      end if
    end do ! end of loop over multiple ridges

    ! Calculate energy change for output to CAM's energy checker.
    call energy_change(dt, p, u, v, ptend%u(:ncol,:), ptend%v(:ncol,:), ptend%s(:ncol,:), de)
    flx_heat(:ncol) = de

    call outfld('TAUARDG'//trim(type)//'X', taurx,  ncol, lchnk)
    call outfld('TAUARDG'//trim(type)//'Y', taury,  ncol, lchnk)

    if (trim(type) == 'BETA') then
      fname(1) = 'TAUGWX'
      fname(2) = 'TAUGWY'
      fname(3) = 'UTGWORO'
      fname(4) = 'VTGWORO'
    else if (trim(type) == 'GAMMA') then
      fname(1) = 'TAURDGGMX'
      fname(2) = 'TAURDGGMY'
      fname(3) = 'UTRDGGM'
      fname(4) = 'VTRDGGM'
    else
      call endrun('gw_rdg_calc: FATAL: type must be either BETA or GAMMA'&
                  //' type= '//type)
    end if

    call outfld(fname(1), taurx(:,pver+1), ncol, lchnk)
    call outfld(fname(2), taury(:,pver+1), ncol, lchnk)
    call outfld(fname(3), utrdg,  ncol, lchnk)
    call outfld(fname(4), vtrdg,  ncol, lchnk)
    call outfld('TTGWORO', ttrdg / cpair,  ncol, lchnk)

    deallocate(tau, gwut, c)

  end subroutine gw_rdg_calc

  ! Add all history fields for a gravity wave spectrum source.
  subroutine gw_spec_addflds(prefix, scheme, band, history_defaults)
  
    use cam_history, only: addfld, add_default, register_vector_field

    ! One character prefix prepended to output fields.
    character(1), intent(in) :: prefix
    ! Gravity wave scheme name prepended to output field descriptions.
    character(*), intent(in) :: scheme
    ! Wave speeds.
    type(GWBand), intent(in) :: band
    ! Whether or not to call add_default for fields output by WACCM.
    logical, intent(in) :: history_defaults

    integer l
    ! 7 chars is enough for "-100.00"
    character(7) fnum
    ! 10 chars is enough for "BTAUXSn32"
    character(10) dumc1x, dumc1y
    ! Allow 80 chars for description
    character(80) dumc2

    ! Overall wind tendencies.
    call addfld (trim(prefix) // 'UTGWSPEC', ['lev'], 'A', 'm/s2', trim(scheme) // ' U tendency - gravity wave spectrum')
    call addfld (trim(prefix) // 'VTGWSPEC', ['lev'], 'A', 'm/s2', trim(scheme) // ' V tendency - gravity wave spectrum')
    call register_vector_field(trim(prefix) // 'UTGWSPEC', trim(prefix) // 'VTGWSPEC')

    call addfld (trim(prefix) // 'TTGWSPEC', ['lev'], 'A', 'K/s' , trim(scheme) // ' T tendency - gravity wave spectrum')

    ! Wind tendencies broken across five spectral bins.
    call addfld(trim(prefix) // 'UTEND1', ['lev'], 'A','m/s2', trim(scheme) // ' U tendency   c < -40')
    call addfld(trim(prefix) // 'UTEND2', ['lev'], 'A','m/s2', trim(scheme) // ' U tendency  -40 < c < -15')
    call addfld(trim(prefix) // 'UTEND3', ['lev'], 'A','m/s2', trim(scheme) // ' U tendency  -15 < c <  15')
    call addfld(trim(prefix) // 'UTEND4', ['lev'], 'A','m/s2', trim(scheme) // ' U tendency   15 < c <  40')
    call addfld(trim(prefix) // 'UTEND5', ['lev'], 'A','m/s2', trim(scheme) // ' U tendency   40 < c ')

    ! Reynold's stress toward each cardinal direction, and net zonal stress.
    call addfld(trim(prefix) // 'TAUE'  , ['ilev'], 'A','Pa', trim(scheme) // ' Eastward Reynolds stress')
    call addfld(trim(prefix) // 'TAUW'  , ['ilev'], 'A','Pa', trim(scheme) // ' Westward Reynolds stress')
    call addfld(trim(prefix) // 'TAUNET', ['ilev'], 'A','Pa', trim(scheme) // ' E+W Reynolds stress')
    call addfld(trim(prefix) // 'TAUN'  , ['ilev'], 'A','Pa', trim(scheme) // ' Northward Reynolds stress')
    call addfld(trim(prefix) // 'TAUS'  , ['ilev'], 'A','Pa', trim(scheme) // ' Southward Reynolds stress')

    ! Momentum flux in each direction.
    call addfld(trim(prefix) // 'EMF', ['lev'], 'A','Pa', trim(scheme) // ' Eastward MF')
    call addfld(trim(prefix) // 'WMF', ['lev'], 'A','Pa', trim(scheme) // ' Westward MF')
    call addfld(trim(prefix) // 'NMF', ['lev'], 'A','Pa', trim(scheme) // ' Northward MF')
    call addfld(trim(prefix) // 'SMF', ['lev'], 'A','Pa', trim(scheme) // ' Southward MF')

    ! Temperature tendency terms.
    call addfld(trim(prefix) // 'TTGWSDF', ['lev'], 'A', 'K/s', &
      trim(scheme)//' t tendency - diffusion term')
    call addfld(trim(prefix) // 'TTGWSKE', ['lev'], 'A', 'K/s', &
      trim(scheme) // ' t tendency - kinetic energy conversion term')

    ! Gravity wave source spectra by wave number.
    do l = -band%ngwv, band%ngwv
      ! String containing reference speed.
      write(fnum, '(f7.2)') band%cref(l)

      dumc1x = tau_fld_name(l, prefix, x_not_y=.true.)
      dumc1y = tau_fld_name(l, prefix, x_not_y=.false.)
      dumc2 = trim(scheme) // ' tau at c= ' // trim(fnum) // ' m/s'
      call addfld(trim(dumc1x), ['lev'], 'A', 'Pa', dumc2)
      call addfld(trim(dumc1y), ['lev'], 'A', 'Pa', dumc2)
    end do

    if (history_defaults) then
      call add_default(trim(prefix) // 'UTGWSPEC', 1, ' ')
      call add_default(trim(prefix) // 'VTGWSPEC', 1, ' ')
      call add_default(trim(prefix) // 'TTGWSPEC', 1, ' ')
      call add_default(trim(prefix) // 'TAUE'    , 1, ' ')
      call add_default(trim(prefix) // 'TAUW'    , 1, ' ')
      call add_default(trim(prefix) // 'TAUNET'  , 1, ' ')
      call add_default(trim(prefix) // 'TAUN'    , 1, ' ')
      call add_default(trim(prefix) // 'TAUS'    , 1, ' ')
    end if

  end subroutine gw_spec_addflds

  ! Outputs for spectral waves.
  subroutine gw_spec_outflds(prefix, lchnk, ncol, band, c, u, v, xv, yv, &
                             gwut, dttdf, dttke, tau, utgw, vtgw, ttgw, taucd)

    use gw_common, only: west, east, south, north

    ! One-character prefix prepended to output fields.
    character(1), intent(in) :: prefix
    ! Chunk and number of columns in the chunk.
    integer, intent(in) :: lchnk
    integer, intent(in) :: ncol
    ! Wave speeds.
    type(GWBand), intent(in) :: band
    ! Wave phase speeds for each column.
    real(r8), intent(in) :: c(ncol,-band%ngwv:band%ngwv)
    ! Winds at cell midpoints.
    real(r8), intent(in) :: u(ncol,pver)
    real(r8), intent(in) :: v(ncol,pver)
    ! Unit vector in the direction of wind at source level.
    real(r8), intent(in) :: xv(ncol)
    real(r8), intent(in) :: yv(ncol)
    ! Wind tendency for each wave.
    real(r8), intent(in) :: gwut(ncol,pver,-band%ngwv:band%ngwv)
    ! Temperature tendencies from diffusion and kinetic energy.
    real(r8) :: dttdf(ncol,pver)
    real(r8) :: dttke(ncol,pver)
    ! Wave Reynolds stress.
    real(r8), intent(in) :: tau(ncol,-band%ngwv:band%ngwv,pver)
    ! Zonal and meridional total wind tendencies.
    real(r8), intent(in) :: utgw(ncol,pver)
    real(r8), intent(in) :: vtgw(ncol,pver)
    ! Temperature tendencies.
    real(r8), intent(in) :: ttgw(ncol,pver)
    ! Reynolds stress for waves propagating in each cardinal direction.
    real(r8), intent(in) :: taucd(ncol,pver+1,4)

    ! Indices
    integer i, k, l
    integer ix(ncol, -band%ngwv:band%ngwv), iy(ncol,-band%ngwv:band%ngwv)
    integer iu(ncol), iv(ncol)

    ! Zonal wind tendency, broken up into five bins.
    real(r8) utb(ncol,pver,5)
    ! Definition of the bin boundaries.
    real(r8), parameter :: bounds(4) = [ -40.0_r8, -15.0_r8, 15.0_r8, 40.0_r8 ]

    ! Momentum flux in the four cardinal directions.
    real(r8) mf(ncol,pver,4)

    ! Wave stress in zonal/meridional direction
    real(r8) taux(ncol,-band%ngwv:band%ngwv,pver)
    real(r8) tauy(ncol,-band%ngwv:band%ngwv,pver)

    ! Temporaries for output
    real(r8) dummyx(ncol,pver)
    real(r8) dummyy(ncol,pver)
    ! Variable names
    character(10) dumc1x, dumc1y

    ! Accumulate wind tendencies binned according to phase speed.

    utb = 0

    ! Find which output bin the phase speed corresponds to.
    ix = find_bin(c)

    ! Put the wind tendency in that bin.
    do l = -band%ngwv, band%ngwv
      do k = 1, pver
        do i = 1, ncol
          utb(i,k,ix(i,l)) = utb(i,k,ix(i,l)) + gwut(i,k,l)
        end do
      end do
    end do

    ! Find just the zonal part.
    do l = 1, 5
      do k = 1, pver
        utb(:,k,l) = utb(:,k,l) * xv
      end do
    end do

    call outfld(trim(prefix) // 'UTEND1', utb(:,:,1), ncol, lchnk)
    call outfld(trim(prefix) // 'UTEND2', utb(:,:,2), ncol, lchnk)
    call outfld(trim(prefix) // 'UTEND3', utb(:,:,3), ncol, lchnk)
    call outfld(trim(prefix) // 'UTEND4', utb(:,:,4), ncol, lchnk)
    call outfld(trim(prefix) // 'UTEND5', utb(:,:,5), ncol, lchnk)

    ! Output temperature tendencies due to diffusion and from kinetic energy.
    call outfld(trim(prefix) // 'TTGWSDF', dttdf / cpair, ncol, lchnk)
    call outfld(trim(prefix) // 'TTGWSKE', dttke / cpair, ncol, lchnk)

    ! Output tau broken down into zonal and meridional components.
    taux = 0
    tauy = 0
    ! Project c, and convert each component to a wavenumber index.
    ! These are mappings from the wavenumber index of tau to those of taux
    ! and tauy, respectively.
    do l = -band%ngwv, band%ngwv
      ix(:,l) = c_to_l(c(:,l) * xv)
      iy(:,l) = c_to_l(c(:,l) * yv)
    end do
    ! Find projection of tau.
    do k = 1, pver
      do l = -band%ngwv, band%ngwv
        do i = 1, ncol
          taux(i,ix(i,l),k) = taux(i,ix(i,l),k) + abs(tau(i,l,k) * xv(i))
          tauy(i,iy(i,l),k) = tauy(i,iy(i,l),k) + abs(tau(i,l,k) * yv(i))
        end do
      end do
    end do

    do l = -band%ngwv, band%ngwv
      dummyx = taux(:,l,:)
      dummyy = tauy(:,l,:)

      dumc1x = tau_fld_name(l, prefix, x_not_y=.true.)
      dumc1y = tau_fld_name(l, prefix, x_not_y=.false.)

      call outfld(dumc1x, dummyx, ncol, lchnk)
      call outfld(dumc1y, dummyy, ncol, lchnk)
    end do

    ! Output momentum flux in each cardinal direction.
    mf = 0
    do k = 1, pver
      ! Convert wind speed components to wavenumber indices.
      iu = c_to_l(u(:,k))
      iv = c_to_l(v(:,k))

      ! Sum tau components in each cardinal direction.
      ! Split west/east and north/south based on whether wave speed exceeds
      ! wind speed.
      do l = -band%ngwv, band%ngwv
        where (iu > l)
          mf(:,k,west) = mf(:,k,west) + taux(:,l,k)
        else where
          mf(:,k,east) = mf(:,k,east) + taux(:,l,k)
        end where

        where (iv > l)
          mf(:,k,south) = mf(:,k,south) + tauy(:,l,k)
        else where
          mf(:,k,north) = mf(:,k,north) + tauy(:,l,k)
        end where
      end do
    end do

    call outfld(trim(prefix) // 'WMF', mf(:,:,west ), ncol, lchnk)
    call outfld(trim(prefix) // 'EMF', mf(:,:,east ), ncol, lchnk)
    call outfld(trim(prefix) // 'SMF', mf(:,:,south), ncol, lchnk)
    call outfld(trim(prefix) // 'NMF', mf(:,:,north), ncol, lchnk)

    ! Simple output fields written to history file.
    ! Total wind tendencies.
    call outfld(trim(prefix) // 'UTGWSPEC', utgw, ncol, lchnk)
    call outfld(trim(prefix) // 'VTGWSPEC', vtgw, ncol, lchnk)
    call outfld(trim(prefix) // 'TTGWSPEC', ttgw, ncol, lchnk)
    ! Tau in each direction.
    call outfld(trim(prefix) // 'TAUE'    , taucd(:,:,east ), ncol, lchnk)
    call outfld(trim(prefix) // 'TAUW'    , taucd(:,:,west ), ncol, lchnk)
    call outfld(trim(prefix) // 'TAUN'    , taucd(:,:,north), ncol, lchnk)
    call outfld(trim(prefix) // 'TAUS'    , taucd(:,:,south), ncol, lchnk)
    call outfld(trim(prefix) // 'TAUNET'  , taucd(:,:,east)+taucd(:,:,west), ncol, lchnk)

  contains

    ! Given a value, finds which bin marked by "bounds" the value falls into.
    elemental function find_bin(val) result(idx)

      real(r8), intent(in) :: val

      integer idx

      ! We just have to count how many bounds are exceeded.
      if (val >= 0) then
        idx = count(val > bounds) + 1
      else
        idx = count(val >= bounds) + 1
      end if

    end function find_bin

    ! Convert a speed to a wavenumber between -ngwv and ngwv.
    elemental function c_to_l(c) result(l)

      real(r8), intent(in) :: c

      integer l

      l = min( max(int(c/band%dc),-band%ngwv), band%ngwv )

    end function c_to_l

  end subroutine gw_spec_outflds

  ! Generates names for tau output across the wave spectrum (e.g.
  ! BTAUXSn01 or TAUYSp05).
  ! Probably this should use a wavenumber dimension on one field rather
  ! than creating a ton of numbered fields.
  character(9) pure function tau_fld_name(l, prefix, x_not_y)

    integer, intent(in) :: l           ! Wavenumber
    character(1), intent(in) :: prefix ! Single-character prefix for output
    logical, intent(in) :: x_not_y

    character(2) num_str

    tau_fld_name = trim(prefix)

    tau_fld_name = trim(tau_fld_name) // 'TAU'

    if (x_not_y) then
      tau_fld_name = trim(tau_fld_name) // 'XS'
    else
      tau_fld_name = trim(tau_fld_name) // 'YS'
    end if

    if (l < 0) then
      tau_fld_name = trim(tau_fld_name) // 'n'
    else
      tau_fld_name = trim(tau_fld_name) // 'p'
    end if

    write(num_str, '(I2.2)') abs(l)

    tau_fld_name = trim(tau_fld_name) // num_str

  end function tau_fld_name

end module gw_drag
